# I load the needed libraries
library(dplyr)
library(scales)
library(GoFKernel)
library(mvtnorm)
library(gplots)
options(warn=-1)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
Loading required package: KernSmooth
KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009
Attaching package: ‘gplots’
The following object is masked from ‘package:stats’:
lowess
I load the functions from the class file:
source("class_MCMC.R")
I define then the function that I want to use as output of the MCMCs:
# Function to sampled from: n-dim gaussian with chosen sigmas and centers
# posterior_g_inhom = function (theta) {
# sigmas = c(1:length(theta))
# centers = c(seq(length(theta), 1))
# product = 1
# for (i in 1:length(theta)) {
# product = product * exp(-(theta[i] - centers[i])**2/sigmas[i]**2)
# }
# return (product)
# }
cauchy1_gauss2 = function (theta) {
sigmas = c(2.5, 4.3)
centers = c(0.4, 9)
product = 1
for (i in 1:2) {
product = product * exp(-(theta[i] - centers[i])**2/sigmas[i]**2)
}
product = product * (dcauchy(theta[3], -10, 2) + 4*dcauchy(theta[3], 10, 4))
return (product)
}
chosen_function = cauchy1_gauss2
Then I only have to determine the parameters for the initialization = the "hyperparameters" of the simulations
# The initial parameters are:
init = c(2, 4, 10)
std = diag(1, 3)
N = as.integer(1e5)
burn_in = as.integer(1e4)
print_step = as.integer(1e2)
# print_init = as.integer(1e3)
N_tot = N + burn_in
# For Haario:
epsilon = 0.001
# MVTNORM
# Evaluate then the MCMC
mcmc_g = random_steps_mvtnorm (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE)
# Selecting the sequence after the burn-in
mcmc_g = mcmc_g[burn_in:N, ]
# Plotting the results
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 76.43091 %
# MVTNORM GIBBS
mcmc_g = random_steps_mvtnorm_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 83.94091 %
# # SIMPLE ADAPTIVE
# mcmc_g = random_steps_simple (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
# gamma_function = gamma_series_exp, halved_step = burn_in)
# mcmc_g = mcmc_g[burn_in:N, ]
# show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
# # SIMPLE ADAPTIVE GIBBS
# mcmc_g = random_steps_simple_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
# gamma_function = gamma_series_exp, halved_step = burn_in)
# mcmc_g = mcmc_g[burn_in:N, ]
# show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
# HAARIO
mcmc_g = random_steps_haario (func_wanted = chosen_function, theta_init = init, n_samples = N_tot,
sigma = std, print_accept=TRUE, t_0 = burn_in, eps = epsilon)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 16.50455 %
Final mean = 0.4193031 9.023863 5.946447
Final covariance matrix =
[,1] [,2] [,3]
[1,] 6.125695 11.89185 7.440331
[2,] 11.891855 295.06024 162.545150
[3,] 7.440331 162.54515 894.842882
# HAARIO GIBBS
mcmc_g = random_steps_haario_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot,
sigma = std, print_accept=TRUE, t_0 = burn_in, eps = epsilon)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 34.34333 %
Final mean = 0.3690438 8.992309 6.83135
Final covariance matrix =
[,1] [,2] [,3]
[1,] 6.470395 11.43808 8.612761
[2,] 11.438083 293.97645 230.063636
[3,] 8.612761 230.06364 830.000087
# RAO
mcmc_g = random_steps_AM_rao (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
gamma_function = gamma_series_exp, halved_step = burn_in/2)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 43.51091 %
Final mean = 0.3750626 9.02413 6.453105
Final covariance matrix =
[,1] [,2] [,3]
[1,] 3.13846802 -0.04642791 0.222659
[2,] -0.04642791 9.32302506 -1.267430
[3,] 0.22265900 -1.26743030 405.974801
# RAO GIBBS
mcmc_g = random_steps_AM_rao_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
gamma_function = gamma_series_exp, halved_step = burn_in/2)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 59.94212 %
Final mean = 0.4077146 9.051977 6.745132
Final covariance matrix =
[,1] [,2] [,3]
[1,] 3.690542211 0.002453223 -0.4264274
[2,] 0.002453223 8.570809601 0.5417770
[3,] -0.426427380 0.541777006 269.2037717
# GLOBAL
mcmc_g = random_steps_global (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
gamma_function = gamma_series_exp, halved_step = burn_in)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 55.27545 %
Final mean = 0.556769 9.134921 5.749301
Final lambda = -0.5730512
Final covariance matrix =
[,1] [,2] [,3]
[1,] 2.8205581 -0.2656854 -1.682357
[2,] -0.2656854 8.9788965 1.576610
[3,] -1.6823573 1.5766099 219.212919
# GLOBAL GIBBS
mcmc_g = random_steps_global_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
gamma_function = gamma_series_exp, halved_step = burn_in)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 69.36121 %
Final mean = 0.1872337 9.060671 6.700329
Final lambda = -0.9273487
Final covariance matrix =
[,1] [,2] [,3]
[1,] 3.52341291 -0.1182699 -0.02448046
[2,] -0.11826987 8.3311801 0.98728850
[3,] -0.02448046 0.9872885 138.62285448
